%导言区
\documentclass{article}%文件格式
\usepackage{ctex,amsmath}%默认引入的包
\usepackage{url}%bibtex参考网页需要的包
\usepackage{graphicx}%引入图片需要的包
\usepackage{float}%修改图片大小scale需要的包
\usepackage{hyperref}%使参考资料能够变成超链接需要的包
\usepackage{listings}%插入代码格式化需要的包
\usepackage{color}%修改文字颜色需要的包
\usepackage[framed,numbered,autolinebreaks,useliterate]{mcode}%美观
\usepackage{cite}%引用参考文献bitex需要的包
\def\degree{${}^{\circ}$}%度数的角标
\graphicspath{{figures/}}%设置图片引用路劲为当前文件夹下面的figures文件夹

\numberwithin{equation}{section}%设置公式后面的编号加上章节
\numberwithin{figure}{section}%设置图片后面的编号加上章节
\hypersetup{
	colorlinks=true,
	linkcolor=blue,
	filecolor=blue,      
	urlcolor=blue,
	citecolor=cyan,
}%设置超链接格式

%% 设置文件整体版面布局
\addtolength{\topmargin}{-54pt}
\setlength{\oddsidemargin}{0.63cm}  % 3.17cm - 1 inch
\setlength{\evensidemargin}{\oddsidemargin}
\setlength{\textwidth}{14.66cm}
\setlength{\textheight}{24.00cm}    % 24.62

\title{QPSK收发机工作总结}%文件标题
\date{}


%正文区
\begin{document}
\maketitle
\tableofcontents
\newpage
\section{项目简介}


\section{需求分析}
\subsection{技术需求}
\begin{enumerate}
	\item 符号速率：2.5MHz
	\item 发射机DAC时钟频率：60MHz
\end{enumerate}

\subsection{功能需求}
\begin{enumerate}
	\item {\color{red} 缺}
	\item {\color{red} 缺}
\end{enumerate}

\section{总体设计}
{\color{red} 参考\cite{丁鹏仁2015基于,黄凌2010基于}确定发射机和接收机的整体架构，这两篇文献讲的很全}

QPSK收发机整体上可以分为QPSK发射机、信道仿真、QPSK接收机3个部分，接下来将分别介绍3个部分的构成以及功能。

QPSK发射机结构如下图所示，整体上分为3个部分，分别是将比特流映射为IQ通道符号流的符号映射、对符号流限制带外泄露以及去除码间串扰形成基带信号的成型滤波、对基带信号调频到中频上的正交调制。
\begin{figure}[H]
	\centering
	\includegraphics[scale=0.5]{modulation.png}
	\caption{QPSK发射机结构图}
	\label{fig:sender structure}
\end{figure}

发射机通过天线发射出去的信号，经过信道传输，最后通过接收机的天线接收下来进行信号处理。接收机总体上分为4个部分，首先是正交解调将信号从中频搬移到基频上，然后下采样将采样速率从60MHz降低到10MHz，匹配滤波以增大信噪比，相干解调解调出基带信号，均衡消除非理想信道传输的影响，最后对符号进行判决，确定最后接收到的信息。其结构如下所示：
\begin{figure}[H]
	\centering
	\includegraphics[scale=0.6]{receiver structure.png}
	\caption{QPSK接收机结构图}
\end{figure}

\section{软件设计}
\subsection{QPSK介绍}\label{chap1:qpsk introduction}
在介绍QPSK之前，我们先来介绍一下为什么需要QPSK带通信号调制方法。这部分我们将探讨基带调制，然后在此基础上讨论带通信号调制，讨论带通信号调制各种方法的区别以及他们的适用性。
\subsubsection{基带信号传输}\label{chap1.1: baseband modulation}
首先我们先介绍基带信号传输，以一个最简单的例子为例，假设我们现在期望将"Hello World!"这一串字符发送到数公里以外的某台设备上，熟悉二进制计算机系统的可能会很快想到我们可以将字符串转换为01序列利用电平信号的高低传输信息。这就是二进制基带信号传输的朴素道理，其中将01序列转换为电平表示这一过程就称为\textbf{脉冲调制。}

如下图所示，我们令每个数字对应的电平持续时间为$T$，则我们可以如(b)所示将1表示成为持续时间内有脉冲，将0表示为持续时间内没有脉冲。如果我们将脉冲的时间增大到与持续时间$T$相等，那么就如(c)所示，此时我们将1表示成为持续时间为$T$的高电平，将0表示为持续时间为$T$的低电平。这种脉冲调制得到的二进制电平波形被称为脉冲编码调制(Pulse-Code Modulation, PCM)波形。\textbf{如下图(c)所示的就是最为常见的不归零(NRZ)波形}\cite{sklar2015数字通信}。
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.6]{chap1_pcm.png}
    \caption{表示二进制数字信息的波形示例}
    \label{fig:chap1_pcm}
\end{figure}

\subsubsection{带通信号传输与调制}\label{chap1.2:bandpass modulation}
数字调制是将数字符号转换为适合信道特性的波形的过程。在\ref{chap1.1: baseband modulation}中我们介绍了基带信号调制，可以看到它的波形通常具有整形脉冲的形式，而在带通调制(bandpass Modulation)中则是利用整形脉冲去调制正弦信号，这个正弦信号称为载波(carrier)。\textbf{将载波转换为电磁场就可以实现无线通信。}

\textbf{要实现基带信号的无线传输必须通过载波}。因为电磁场的传播必须利用天线，而天线的尺寸和信号的波长成正比(一般天线的长度为波长的四分之一)，因此对于一个频率为3000Hz的基带信号而言，需要的天线尺寸为$2.5\times10^8m$，这显然是不现实的，因此我们可以将该基带信号调制到频率为900MHz的载波上，此时等效天线尺寸则为$8cm$。

因此，对于各种无线传输系统来看，利用载波进行带通调制是非常有必要的。除了可以降低天线尺寸以外，还可以实现频分复用、扩频调制等实现更为复杂有效的调制方法。

\subsubsection{QPSK正交调制}\label{chap1.3:qpsk modulation}
正交相移键控(Quadrature Phase Shift Keying, QPSK)是目前最为常用的一种卫星数字调制方式，它具有较高的频率利用率、较强的抗干扰性以及容易实现\cite{黄凌2010基于}。对于所有的相移键控(Phase Shift Keying, PSK)调制方法，都可以用下式表示
\begin{equation}\label{eq:qpsk modulation}
        s(t) = \sum_n a_ng(t-nT_s)\cos(w_ct+\phi_n)
\end{equation}
其中$a_n$为二进制数字，对于QPSK，$a_n$的取值一共有4个。$T_s$为每个码元的持续时间，$w_c$为载波角频率，$\phi_n$为第$n$个码元的载波相位取值。也就是说我们使用码元改变当前载波在持续时间$T_s$内的相位。如下图所示：
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.9]{chap1.3_qpsk_waveform.jpg}
    \caption{QPSK调制波形}
    \label{fig:chap1.3_qpsk_waveform}
\end{figure}

相比于二进制相移键控调制(Binary Phase Shift Keying,BPSK)只利用了0\degree 和180\degree 两个载波相位，QPSK利用了4个载波相位。那么很自然的，我们想到可以修改载波的相位生成发射信号，这是QPSK常用的相位选择法。但是我们这里将从复包络入手介绍实信号的复数表示，从而引出调相法。

对于任意一个\textbf{实带通信号s(t)}，我们可以将其表示为下面的复数形式：
\begin{equation}
    s(t) = \Re\{g(t)e^{jw_ct}\}
\end{equation}
其中$g(t)$为复包络，其表达式为
\begin{equation}
    g(t) = x(t) + jy(t) = |g(T)|e^{j\theta(t)} = R(t)e^{j\theta(t)}
\end{equation}
其中$g(t)$是复数形式的基带信息或数据，$e^{jw_ct}$是载波的复数表示。将这两项相乘就是调制，\textbf{乘积的实部s(t)就是发射信号}，所以有
\begin{equation}
    \begin{split}
        s(t) &= \Re\{g(t)e^{jw_ct}\}\\
        &= \Re\{[x(t)+jy(t)][\cos w_ct+j\sin w_ct]\}\\
        &= x(t)\cos w_ct - y(t)\sin w_ct
    \end{split}
\end{equation}
因此对于基带信号，我们将其分解(串/并变换)为同相信号$x(t)$和正交信号$y(t)$，然后分别和余弦信号以及载波信号相乘，最后作差即可得到发射信号。其结构如下图所示：
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.3]{chap1.3_qpsk_modulation_structure.png}
    \caption{QPSK正交调制结构示意图}
    \label{fig:chap1.3_qpsk_modulation_structure}
\end{figure}

\subsection{成型滤波}\label{chap2:shaping filter}
在\ref{chap1.3:qpsk modulation}中我们将基带信号分解进行QPSK调制到载波上作为发送信号，但是我们需要注意到基带信号的频谱是很大的，在频谱资源有限的工作条件下可能我们必须对基带信号的频谱进行一定的约束。考虑到基带信号由无数个持续时间为$T_s$的方波组成，我们将其表示为
\begin{equation}
    s(t) = \sum_n a_ng(t-nT_s)
\end{equation}
其中$a_n$为需要发送的二进制数字，考虑到上式是一个随机过程，我们并不知道当前发送的二进制数字是什么，因此我们求其能量谱密度(Power Spectrum Density, PSD)有：
\begin{equation}
    S_{QPSK}(f) = sinc^2(2T_s(f-f_c))+sinc^2(2T_s(f+f_c))
\end{equation}
相比于BPSK的能量谱密度PSD
\begin{equation}
    S_{BPSK}(f) = sinc^2(T_s(f-f_c))+sinc^2(T_s(f+f_c))
\end{equation}
可以看出QPSK的主瓣带宽只有BPSK主瓣带宽的一半，如下图所示。但是及时这样，我们同样可以观察到PSD的旁瓣同样集中了很大的能量，占用了很大的带宽资源，因此我们需要\textbf{设计一个低通滤波器来降低旁瓣的能量，减小带宽的占用}。
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.7]{chap2.1_qpsk_psd.png}
    \caption{QPSK和BPSK的能量谱密度对比示意图}
    \label{fig:chap2.1_qpsk_psd}
\end{figure}

同时考虑到如下图所示的典型数字通信系统的滤波问题。在整个通信系统中存在着各种各样的滤波器，他们对信号造成了不同的影响。而由于系统的滤波作用，接收脉冲之间会发生交迭。脉冲出现拖尾占据了相邻码元间隔，从而干扰了信号检测过程，进而造成误差性能的降低，这类干扰被称为码间串扰(Intersymbol Interference,ISI)。即便此时没有噪声，依旧会发生ISI。
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.4]{chap2.1_isi.png}
    \caption{检测过程中的码间串扰}
    \label{fig:chap2.1_isi}
\end{figure}

综上，我们希望设计一个滤波器，不仅能够起到低通滤波器的作用限制信号旁瓣的能量降低频谱资源占用，还能降低码间串扰对信号的影响。\textbf{奈奎斯特证明了要使码元速率为$R_s$码元/s的信号不存在码间串扰，理论上所需的最小的系统带宽为$R_s/2Hz$，当且仅当系统传输函数$H(f)$如下图左示的矩形函数。}对于基带系统，$H(f)$是单边带宽为$1/2T_s$的矩形函数，系统的冲击响应即为$H(f)$的反傅里叶变换$h(t)=sinc(t/T)$，如下图右所示。

可以观察到，相邻的两个间隔为$T_s$的峰恰好位于彼此的旁瓣处，只要采样正确就不会存在码间串扰。同时奈奎斯特还证明了对于基带系统，无码间串扰地监测马元速率为$1/T$的脉冲所需要的带宽为$1/2T$；换言之，带宽$W=1/2T=R_s/2Hz$的系统在保证无码间串扰的条件下能够支持的最大传输速率为$2W=1/T=R_s$码元/s。对于奈奎斯特脉冲可以表示为$sinc(t/T)$函数和另一个时间函数的乘积，而常见的成型滤波器为升余弦成型滤波以及根升余弦成型滤波器\cite{sklar2015数字通信}。
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.6]{chap2.1_no_isi.png}
    \caption{无码间串扰的奈奎斯特信道}
    \label{fig:chap2.1_no_isi}
\end{figure}

符号映射出来的结果$I_n$依旧是离散的，反映在时域上就是冲激函数序列，显然其频谱是无限宽的。而由于每一个无线通信系统都工作在一定的频带之内，因此必须要对符号映射出来的结果进行约束限制其带外泄露。假设成型滤波器表示为$g(t)$，那么对成型滤波器$g(t)$提出的第一个要求就是它应该是一个低通滤波器。

除此以外，假设信道的响应为$\delta (t)$，接收到的基带信号就是发射时的基带信号，那么有$$r(t)=s(t)=\sum_n I_n g(t-n T_s)$$如果接收机在$t=nT_s$时刻采样，那么希望得到的采样值就应该是符号$I_n$，因此要求$g(t)$有如下的性质：
\begin{equation}
	g(n T_s)=\begin{cases}
		1, & n = 0\\
		0, & \text{else}
	\end{cases}
\end{equation}
也就是说，经过成型之后的波形在$T_s$时的采样点不能被之前的符号影响，这个叫做无码间串扰要求。

综上，成型滤波器有两点要求，首先要求其是一个低通滤波器，其次要求它具有无码间串扰特性。很容易想到的是理想低通滤波器，它的时域响应为$sinc()$函数，很好地满足了无码间串扰的要求，但是它的时域响应衰减太慢了，一但采样时刻稍微偏差，符号的码间串扰就会累加起来造成很大的误差。因此这里将介绍衰减相对理想低通滤波器快速很多的常用的升余弦滤波器(raised cosine filter)。其公式为：
\begin{equation}
	g(t)=\frac{\cos(\pi\alpha\frac{t}{T_s})}{1-(2\alpha \frac{t}{T_s})^2} \frac{\sin \pi \frac{t}{T_s}}{\pi \frac{t}{T_s}}
\end{equation}
可以从公式中看到，右边依然是一个$sinc()$函数的样子，只不过在左边加入了一个分母为$\frac{1}{T_s^2}$的衰减因子，一般来讲3个周期左右时域响应就会衰减到0左右。它的时域以及频域响应如下图所示：
\begin{figure}[H]
	\centering
	\includegraphics[scale=0.6]{rcos.jpg}
	\caption{升余弦滤波器时域频域响应图}
	\label{fig:rcos}
\end{figure}
可以看到时域上它满足无码间串扰，在除了自身采样点之外的采样点的采样值都为0，大致呈现出$sinc()$函数的样子，同时相比$sinc()$函数具有更小的旁瓣。在频域上整体呈现出一个低通滤波器的性质，截止频率为$\frac{1+\alpha}{2T_b}$，其中$\alpha$为滚降因子，主要控制升余弦滤波器的滚降程度，也就是频域上衰减带的宽度。

但是上图所示的升余弦滤波器是非因果的，因此为了实现因果的升余弦滤波器，必须将时域上的结果加上一定的延时，使时域波形整体右移至正半轴上。同时考虑到不可能模拟无限长的升余弦滤波器，在MATLAB仿真的时候必须截断用有限长度的时域点数来近似表示升余弦滤波器的时域波形。

MATLAB中设计升余弦滤波器的函数为$rcosdesign(rollof,span,sps)$函数，其中rollof就是公式中的滚降因子，取值为0-1，而span是指整个升余弦滤波器覆盖的符号个数，也就是每次处理几个符号，\textbf{sps指的是每个符号对应的采样点数}。通常对符号进行4倍过采样，也就是sps=4，一般rollof取0.3-0.8的值都可以，rollof的值越低升余弦滤波器的频域响应越像一个理想的低通滤波器。而升余弦滤波器的延时可以通过MATLAB的grpdelay函数得出，同时从升余弦滤波器的点数上可以看出，一共需要处理$span*sps+1$个采样点，因此延时为$(span*sps)/2$。

因此，如下图所示，符号串$I_n$经过升余弦滤波器成型之后的基带信号波形可以理解为多个升余弦滤波器时域响应的叠加，并且在每个采样时刻其他信号成型结果刚好为0。
\begin{figure}[H]
	\centering
	\includegraphics[scale=0.9]{rcos shaping.eps}
	\caption{升余弦成型时域图}
	\label{fig:rcos shaping}
\end{figure}
但是实际上我们使用的更多的是根升余弦滤波器(root raised cosine filter)，相比升余弦滤波器做了一个取平方根的操作。这是为了方便后面进行匹配滤波而设计的。根升余弦滤波器的特点在于在其他采样点时它的值不为0，但是与另一个根升余弦滤波器级联之后可以消除这个影响。

{\color{red} Q\&A:如果在不使用根升余弦滤波器而是在接收端和发射端都使用升余弦滤波器会对信号带来什么样的影响？}


\subsection{上采样}\label{subsec:software upsample}
{\color{red} 为什么需要上采样——介绍OSR带来的好处（降低抗混叠滤波器的设计难度，增大信噪比等，参考\cite{sklar2015数字通信}的2.4.3和\cite{丁鹏仁2015基于}）；介绍内插滤波器，主要介绍使用什么样的方法在信号之间完成内插，拉格朗日法可以简单介绍，然后引出Farrow内插器，参考}


\subsubsection{为什么需要上采样}
为了使符号串通过升余弦成型滤波器，我们对每个符号都过采样了4次(sps=4)。考虑到原始的符号速率为4MHz，那么4倍过采样之后的基带信号的速率为16MHz，但是考虑到很多时候硬件DAC的速率并不和信号速率相等，因此有必要对信号进行上采样以适应硬件速率。

\subsubsection{上采样实现——以Farrow内插器为例}
一般常用的上采样就是MATLAB提供的$upsample(data,n)$函数，但是$upsample$函数只能提供整数倍的过采样，对于不满足整数倍关系
的采样频率，$upsample$往往无能为力。哪怕是使用$resample(data,up,down)$先整数倍内插$up$次，然后$down$倍抽取也会存在公约数极大以至于无法实现的问题。比如期望从1MHz的信号速率上采样到1.0001MHz上，就需要先以10001倍内插然后10000倍抽取实现，显然这是一个非常大的工作量，不具有实时性。因此有必要寻找一种新的采样手段来解决这个问题。

上采样的本质是实现速率转化的功能，可以假设他的模型如下所示：
\begin{figure}[H]
	\centering
	\includegraphics[scale=0.5]{interpolation_structure.png}
	\caption{上采样的速率转换模型}
	\label{fig:interpolation structure}
\end{figure}
假设采样时钟周期为$T_s$，符号周期为$T$，其中数字信号$x(mT_s)$经过DAC之后得到一个模拟信号，通过滤波器$h(t)$之后得到一个连续时间的输出
\begin{equation}
	y(t)=\sum_m x(mT_s)h(t-mT_s)
\end{equation}
假设对于$y(t)$在时刻$t=kT_i$重新进行采样，其中$T_i$为内插周期，那么采样之后的结果有：
\begin{equation}\label{eq:interpolation}
	y(kT_i)=\sum_m x(mT_s)h(kT_i-mT_s)
\end{equation}
即使在这个速率转化的模块中出现了DAC，但是只要$x(m)$是已知的，并且知道滤波器$h(t)$的冲击响应以及输入信号的采样时间$T_s$、输出信号的采样时间$T_i$，那么就可以根据上式完全计算出内插值。\textbf{其本质就是，知道内插的时刻，知道怎么利用周围的采样点计算内插值就可以得到内插值。}如下图所示，黑色竖线为已有的原始采样速率下的采样值，而红色虚线则为新的采样速率下的内插出来的结果。
\begin{figure}[H]
	\centering
	\includegraphics[scale=0.7]{interpolation.png}
	\caption{采样点关系图}
\end{figure}
这里给出如下定义，首先定义基本指针$m_k$，它表示当前插值位于第$m_k$个采样值后面
\begin{equation}
	m_k = int[kT_i/T_s]
\end{equation}
其中$int[]$为向下取整运算。然后定义小数间隔，也就是当前内插值和基本指针之间的时间差$\mu_k T_s$
\begin{equation}
	\mu_k = kT_i/T_s - m_k
\end{equation}
也就是说只要知道了$m_k$和$\mu_k$，就知道了内插的具体时刻，则\ref{eq:interpolation}可以写成
\begin{equation}
	y(kT_i)= y((m_k+\mu_k)T_s)=\sum_{i=N_1}^{N_2}x[(m_k-i)T_s]h[(i+\mu_k)T_s]
\end{equation}
上式就是数字内插滤波器的基本方程，$\mu_k$和$m_k$表示了$T_s$和$T_i$之间的调整关系。一般常用的是具有Farrow结构的插值器，这个结构只需要从最接近内插时刻的四个连续输入信号$x[(m_k-1)T_s],x[(m_k)T_s],x[(m_k+1)T_s],x[(m_k+2)T_s]$来计算内插值$y[(m_k+\mu_k)T_s]$。这里也恰好符合我们在升余弦成型滤波中每个符号采样4个点的设置。

\textbf{需要注意的是，在现实中我们不可能提前获取$x[(m_k+1)T_s]$的值，只能用已经获取到的四个过去时刻的连续信号如$x[(m_k-3)T_s],x[(m_k-2)T_s],x[(m_k-1)T_s],x[(m_k)T_s]$来实现，所以实际上内插值为$y[(m_k-2+\mu_k)T_s]$}。也就是说Farrow内插之后的结果相比原本的信号有了2个$T_s$的延时。

Farrow结构的内插器如下所示：
\begin{figure}[H]
	\centering
	\includegraphics[scale=0.6]{farrow.png}
	\caption{Farrow结构内插}
\end{figure}
可以看到上图中一共有3条纵向通路以及一条横向通路，他们的计算公式分别为
\begin{equation}\label{eq:farrow}
\begin{split}
	f_1 &= 0.5x(m)-0.5x(m-1)-0.5x(m-2)+0.5x(m-3)\\
	f_2 &= 1.5x(m-1)-0.5x(m)-0.5x(m-2)-0.5x(m-3)\\
	f_3 &= x(m-2)\\
	y_1(k) &= f_1 u(k) u(k) + f_2 u(k) + f_3\\
\end{split}
\end{equation}

farrow内插结果MATLAB仿真如下图所示：
\begin{figure}[H]
	\centering
	\includegraphics[scale=0.9]{matlab_farrow_comparision.png}
	\caption{farrow内插MATLAB仿真结果}
\end{figure}

\subsection{信道传输}\label{subsec:software multi-path}
%2022.1.14等待重写，内容包括：信道会对信号造成什么样的干扰，常见的awgn信道，rayleigh多径信道
\subsubsection{AWGN信道}
\subsubsection{多径信道RayleighChannel}
对于基带信号$s(t)$，假设有$k$条信道，每条信道对于信号会有强度$\alpha_k$以及延时$\tau_k$的影响，则接收机收到的信号$y(t)$可以表述为$k$条信道的卷积叠加：
\begin{equation}
	\begin{split}
		y(t) &=\sum_{n=1}^k \alpha_n s(t-\tau_n)e^{j2\pi f_c(t-\tau_n)}\\
		&=e^{j2\pi f_c}\sum_{n=1}^k \alpha_n s(t-\tau_n)e^{-j2\pi f_c\tau_n}
	\end{split}
\end{equation}

首先先了解一下MATLAB是怎么实现多径信道的，在MATLAB的communication toolbox中有comm.RayleighChannel函数，可以生成具有Rayleigh分布性质的多径信道，其具体使用方法参考MathWorks文档库。

这里需要注意的是，输入的参数AveragePathGains对应的是瑞利分布的均值，实际上产生的信道增益是符合该均值的瑞利分布的随机值( 复增益)，所以当看到输入的均值增益为0dB而输出的信号发生改变也无需惊讶，可以在参数中输入PathGainsOutputPort=1输出随机的增益大小。常见的使用方法如下代码块所示：
\begin{lstlisting}
%% rayleigh Channel
rayleighChannel = comm.RayleighChannel('SampleRate',fs,...
'PathDelays',pathDelays,...
'AveragePathGains',pathGains,...
'PathGainsOutputPort',1);
[rx,gains] = rayleighChannel(tx);
\end{lstlisting}

\textbf{RayleighChannel的实现方法为：首先根据延时$\tau_k$建立多个冲击函数，与sinc函数卷积，sinc函数的边瓣周期由输入参数的采样频率决定。从连续波形上来看，整个信道响应就是多个不同延时、不同幅度的sinc函数的叠加；然后对该连续波形进行周期为$\frac{1}{SampleRate}$的抽样。}

这里特别需要注意的是，rayleighChannel函数产生的是复增益gains，而rayleighChannel的Impulse Response的y轴坐标为Magnitude，其等同于对gains进行abs()运算取模，如下图紫色线条对应的结果所示。而多径信道叠加起来的信道模型(蓝色点)则是将所有sinc函数叠加结果采样的结果，如果直接将所有的gains求和相加是不正确的。
\begin{figure}[H]
	\centering
	\includegraphics[scale=0.9]{channel impulse response.png}
	\caption{信道增益及抽头系数图}
	\label{fig:channel impulse response}
\end{figure}

\subsection{下采样以及匹配滤波}\label{subsec:software dowmsample and shapping}

{\color{red} 两个目的，第一：下采样是为了什么(减小信号的处理量)，会不会降低信噪比？会不会对后续的信号有影响？第二：匹配滤波的作用是什么？主要参考\cite{2014统计信号处理基础,2009统计信号处理,2006随机信号分析与处理}介绍匹配滤波的作用}

与发射端一样，此时的信号的采样速率为60MHz，这对于后面信号的处理是相当快速的，我们希望在处理之前能够将信号的速率相对降低一点来方便后面系统的设计，这里希望下采样之后的采样速率为10MHz，每个符号对应4个采样点以方便后续匹配滤波。因此需要在输出端使用Farrow结构完成一个抽取，注意这里不同于发射端，这里完成的是抽取的作用。

完成抽取之后需要对信号进行匹配滤波以最大化信号的信噪比，考虑到我们在信号成型的时候使用的是根升余弦滤波器，这里我们也使用根升余弦滤波器以完成信号的匹配。

\subsection{AGC自动增益控制}
受信道、接收机发射机相对距离等因素的影响，进入接收机ADC的信号功率会小于发射机发射的功率，当接收到的信号的功率太低时无法进行后续的计算，因为此时的信号很有可能已经被噪声掩盖了或者信号幅度不足以支撑后续的计算将信号和噪声区分开来。因此我们有必要对信号进行合理地放大，并且为了使接收机能够动态的适应接收到的信号的功率大小，我们期望接收机能够根据接收到的功率大小自动调整接收链路的增益，这个过程就是自动增益控制(Automative Gain Control)过程。

这里介绍最简单的AGC环路：假设接收机接收到的信号在$t$时刻为$r(t)$，那么该时刻信号的能量即为$r(t)r^{*}(t)$，如果我们通过一个长度为$N$的队列保存$t$至$t+N-1$一共$N$个信号的能量，然后求平均作为这一时间段的平均功率$p(N)=\frac{\sum_{t}^{t+N-1}r(t)r^{*}(t)}{N}$，有算法如下：
\begin{itemize}
    \item 如果$p(N)\in[min,max]$，那么增益$K=K$
    \item 如果$p(N)<min$，那么增益$K=K+1$
    \item 如果$p(N)>max$，那么增益$K=K+1$
\end{itemize}
其结构如下图所示
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.9]{AGC-structure.png}
    \caption{AGC结构图}
    \label{fig:AGC structure}
\end{figure}


\subsection{相干解调}
在信号传输的过程中，受到多径信道的干扰以及接收机和发射机之间的相对运动导致的多普勒效应的影响，接收机接收到的信号相比发射机发射的信号会有频偏和相偏，此时接收机接收到的信号可以表示为
\begin{equation}
    r(t)=e^{[j2\pi (f_{c}+f_{dopler}) t + \theta + \phi_t]}
\end{equation}
其中$f_{dopler}$为多普勒效应带来的频偏，$\phi_t$为传输过程中的相偏。如果此时$r(t)$不加上载波恢复消除频偏和相偏进行解调，那么解调后的信号为
\begin{equation}
    r(t)=e^{[j2\pi (f_{dopler}) t + \theta + \phi_t]}
\end{equation}
可以看到信号的相位在多普勒频偏的影响下随着时间的变化会不断发生旋转，在星座图上表示即不是固定的QPSK对应的4个点，而是圆圈状的星座图，如下图所示
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.5]{constellation diagram without synchronization.png}
    \caption{未同步的信号的星座图}
    \label{fig:constellation diagram without synchronization}
\end{figure}

因此，为了消除多普勒效应带来的频偏和相偏(虽然最后依然有可能存在固定相位差)，我们需要对信号进行相干解调，使用锁相环(Phase Lock Loop)对信号进行跟踪，这里将使用BPSK和QPSK常用的Costas环进行载波同步。

同时由于接收机和发射机使用的并不是同一个晶振(事实上只要是不在同一个设备上进行调试基本不可能有一样的晶振)，此时两个晶振之间存在时偏，这个误差使接收端无法正确判断每个码元的起始和截止时刻，导致无法在\textbf{最佳的}时刻进行采样判决以增大信噪比，或者说我们期望在眼图开到最大的地方对信号进行采样作为判决依据。这个时候需要使用位同步解决时偏的影响，此处使用的位同步环路为Gardner位同步环路\cite{黄凌2010基于}。将PLL载波同步和Gardner位同步结合起来就是\textbf{相干解调}，其结构图如下图所示
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.65]{synchronization structure.png}
    \caption{相干解调结构图}
    \label{fig:synchronization structure}
\end{figure}

\subsubsection{Costas环载波同步}
假设此时接收机接收到的经过解调的信号为
\begin{equation}
    r(t)=e^{[j2\pi (f_{d}) t + \theta]}
\end{equation}
我们将其展开为IQ通道的信号表示
\begin{equation}
    \begin{split}
        r(t) &= e^{[j2\pi f_{d} t + \theta]}\\
        &=e^{j\theta}(\cos2\pi f_{d}t + j\sin2\pi f_{d}t)\\
        &=(I(t)+jQ(t))(\cos2\pi f_{d}t + j\sin2\pi f_{d}t)\\
        &=I(t)\cos2\pi f_{d}t - Q(t)\sin2\pi f_{d}t + j[I(t)\sin2\pi f_{d}t + Q(t)\cos2\pi f_{d}t]\\
    \end{split}
\end{equation}
将IQ通道分开表示为
\begin{equation}
    \begin{split}
        I{'}(t) &= I(t)\cos2\pi f_{d}t - Q(t)\sin2\pi f_{d}t\\
        Q{'}(t) &= I(t)\sin2\pi f_{d}t + Q(t)\cos2\pi f_{d}t\\
    \end{split}
\end{equation}
我们期望用Costas环消除误差，也就是令$f_d\Rightarrow0$，此时有
\begin{equation}
    \begin{split}
        I{'}(t) &= I(t) \\
        Q{'}(t) &= Q(t) \\
    \end{split}
\end{equation}

整个Costas环跟常见的锁相环PLL没有什么很大的差别，都由3部分组成，分别为计算误差的\textbf{鉴相器}，调整输出相位控制环路稳定性的\textbf{环路滤波器}以及输出载波\textbf{压控振荡器}。其具体结构如图\ref{fig:synchronization structure}左半部分所示。

需要在这里说明的是不同的鉴相器带来的鉴相结果以及导致的环路最后解调出来的结果是不一致的，这里将分别对几种鉴相器进行介绍，并推导他们的稳定点以及适用的场合。

如图\ref{fig:synchronization structure}所示的鉴相器可以表示为
\begin{equation}
    \begin{split}
        error(t) &= sign(I{'}(t))Q{'}t - sign(Q{'}(t))I{'}(t) \\
    \end{split}
\end{equation}
在符号映射中我们提到有两种映射方式，如果我们选用A映射方式，也就是$I(t)$和$Q(t)$取值都为-1或1，在$f_d\Rightarrow0$时，将上式展开有
\begin{equation}
    \begin{split}
        error(t) &= sign(I{'}(t))Q{'}t - sign(Q{'}(t))I{'}(t) \\
        &= I^2(t)\sin2\pi f_{d}t +I(t)Q(t)\cos2\pi f_{d}t + Q^2(t)\sin2\pi f_{d}t -I(t)Q(t)\cos2\pi f_{d}t  \\
        &= I^2(t)\sin2\pi f_{d}t +   Q^2(t)\sin2\pi f_{d}t \\
        &= A\sin2\pi f_{d}t,A = I^2(t) + Q^2(t) =  \sqrt{2}    \\
    \end{split}
\end{equation}
也就是说选用A映射方式下采用如图\ref{fig:synchronization structure}所示的鉴相器的鉴相结果是一个正弦鉴相器，可以观察到在相位为$\frac{k\pi}{2}$时鉴相器的输出为0达到稳定。\textbf{然而如果使用B映射方法，那么鉴相结果和A映射方法的鉴相结果将会不一样}。在MATLAB中仿真A映射方法和B映射方式下的如图\ref{fig:synchronization structure}所示的鉴相器结果如下图所示
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.7]{pd result with AB mapping.png}
    \caption{不同映射方法对应的鉴相结果示意图}
    \label{fig:pd result with AB mapping}
\end{figure}
可以看到A映射方法在相位为$\frac{k}{2}\pi$处达到稳定，而B映射方法在相位为$\frac{\pi}{4}+\frac{k}{2}\pi$处达到稳定。

除了如图\ref{fig:synchronization structure}所示的鉴相器以外，还有常见的\textbf{点乘}鉴相器，即直接将I通道的信号和Q通道的信号相乘作为鉴相器的输出，可以表示为
\begin{equation}
    \begin{split}
        error(t) &= I{'}(t)Q{'}(t)\\
        &= \frac{[I^2(t)-Q^2(t)]}{2}\sin(4\pi f_d t) + \frac{I(t)Q(t)]}{2}\cos(4\pi f_d t)\\
    \end{split}
\end{equation}
显然这个鉴相器的符号会随着IQ通道的数据发生改变，这显然是不符合鉴相器的需求的，我们在鉴相的时候应该抛开接收到的数据带来的对鉴相器的结果的影响。因此对于QPSK相干解调，为了消除数据符号的影响，使用的是\textbf{平方点乘鉴相器}，即
\begin{equation}
    \begin{split}
        error(t) &= I{'}^2(t)Q{'}^2(t)\\
    \end{split}
\end{equation}
MATLAB仿真结果如下图所示
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.7]{pd result with AB mapping dot.png}
    \caption{不同映射方法对应的\textbf{平方点乘鉴相}结果示意图}
    \label{fig:pd result with AB mapping dot}
\end{figure}

从上面的推导可以看出，在选择不同的映射方法时使用不同的鉴相器会对IQ通道解调出来的数据有不同的影响，是后续进行下一步信息获取处理时需要考虑的问题。

在完成了鉴相器的选择之后，我们需要根据鉴相器的输出性质设计环路滤波器使环路滤波器呈现负反馈的特性，能够不断地减小本地载波信号和实际接收信号之间的相位差。常用的一般为二阶数字环路，其结构如下图所示，由一个积分部分和一个比例部分组成，其本质是一个PI控制器。
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.9]{相干解调-costas-环路滤波器结构.png}
    \caption{环路滤波器结构示意图}
    \label{fig:lf stucture}
\end{figure}
压控振荡器以环路滤波器的输出为输入，改变输出的本地载波的相位。其具体解释为：当鉴相器输出了当前的相位差之后，环路滤波器计算出一个合理的数值控制压控振荡器改变相位步进值，该相位步进值累加到压控振荡器的相位寄存器上作为当前时刻的相位输出。当鉴相器输出为0，即此时无相位差时，环路滤波器在积分环节的作用下会保持一个恒定的值控制压控振荡器输出一个恒定的相位步进值控制当前时刻输出的相位确保鉴相误差始终为0。

因此环路滤波器PI部分的参数设计就非常的重要，此处将不做过多的推导，可以参考\cite{黄凌2010基于,2019全球卫星导航系统原理,2009gps}相关部分内容，这里直接给出计算公式
\begin{equation}
    \begin{split}
        C_1 &= \frac{2\epsilon \omega_n T}{K_d}\\
        C_2 &= \frac{(\omega_n T)^2}{K_d}\\
    \end{split}
\end{equation}
其中$\epsilon$表示环路阻尼系数，在二阶系统中，兼顾超调性和快速性的最优阻尼系数为0.707，而$\omega_n$为环路阻尼震荡频率，$K_d$为环路增益，$T$为进入环路的两个信号之间的时间差。其实普遍来看环路滤波器的参数设计是比较取巧的，可以参考PID参数设计的相关资料。\textbf{一般来看PI的参数给的较小对环路的影响不会很大，只是会影响环路的收敛速度，但是环路终究会收敛的，只要参数在稳定范围内。}

总的来看，整个Costas环载波同步的步骤如下：
\begin{enumerate}
    \item 将环路输入信号和本地载波进行鉴相计算相位差
    \item 使用环路滤波器对相位差进行计算得出相位步进值$\Delta\phi$
    \item 在压控振荡器中令当前输出的信号的相位$\phi$为$\phi_{t-1}+\Delta\phi$
    \item 重复上述过程
\end{enumerate}


\subsubsection{Gardner环位同步}\label{subsec:software gardner}
经过下采样之后的信号采样频率为16MHz，对应每个符号有4个采样点，这里存在一个问题，\textbf{什么时候选取哪一个采样点作为这个符号的样本？}这里将介绍Gardner环位同步方法，该方法有两个特点，首先，Gardner环与载波同步完全独立，可以在载波同步之前也可以在载波同步之后运行，不依赖于载波同步；其次，Gardner环位同步法只要求每个符号2个采样点即可判断出最佳采样点。

Gardner环位同步法期望能够采样在眼图开到最大的位置作为符号的判决依据，对于每个符号4个采样点，Gardner环会通过Farrow内插抽取2个采样点作为这个符号的采样值，这两个点分别为strobe点(出现在符号的峰值时刻)以及midstrobe点(出现在两个符号峰值之间)。

对于QPSK调制，用$y_I(k)$和$y_I(k-1)$表示I通道的两个符号峰值采样点strobe点，用$y_Q(k)$和$y_Q(k-1)$表示Q通道的两个符号峰值采样点strobe点，用$y_I(k-1/2)$表示I通道的两个符号峰值采样点之间的midstrobe点，用$y_Q(k-1/2)$表示Q通道的两个符号峰值采样点之间的midstrobe点。则Gardner环的定时误差检测算法可以表示为
\begin{equation}\label{eq:gardner}
	\mu(k) = y_I(k-1/2)[y_I(k)-y_I(k-1)]+y_Q(k-1/2)[y_Q(k)-y_Q(k-1)]
\end{equation}
如果使用$sign()$函数代替部分运算，则有
\begin{equation}\label{eq:gardner sign}
	\mu(k) = y_I(k-1/2)[sign(y_I(k))-sign(y_I(k-1))]+y_Q(k-1/2)[sign(y_Q(k))-sign(y_Q(k-1))]
\end{equation}
Gardner环通过检测定时误差不断地调整Farrow内插的时间来作为新的strobe和midstrobe点，可以从上式看出，当完美采样的时候，两个符号峰值之间的那个midstrobe点应该为0，那么此时的定时检测误差$\mu$为0，然而当没有完美采样的时候，此时定时检测误差必定不为0，就可以修改Farrow内插的时间更改采样时间。

\textbf{从式\ref{eq:gardner}中可以看出，每次完成了strobe点的采样之后就会计算一次定时误差，但是Gardner环怎么知道当前这个点是midstrobe还是strobe点呢？如果在midstrobe点计算误差会发生什么呢？会不会影响Gardner环的性能呢？}这里将依次列举所有的情况来讨论说明这个问题。
\begin{enumerate}
	\item 在strobe点执行误差检测，相比最佳采样提前
	\begin{figure}[H]
		\centering
		\includegraphics[scale=0.7]{strobe before.png}
		\caption{情况一：在strobe点误差检测，相比完美最佳提前}
	\end{figure}
	如上图所示，易知此时的误差$\mu(k)<0$，也就表明此时相比完美采样提前了，Gardner环位同步会根据这个误差通过负反馈环路增大采样间隔修正下次采样的时间。
	
	\item 在strobe点执行误差检测，相比最佳采样滞后
	\begin{figure}[H]
		\centering
		\includegraphics[scale=0.7]{strobe later.png}
		\caption{情况二：在strobe点误差检测，相比完美最佳滞后}
	\end{figure}
	如上图所示，易知此时的误差$\mu(k)>0$，也就表明此时相比完美采样滞后了，Gardner环位同步会根据这个误差通过负反馈环路减小采样间隔修正下次采样的时间。
	
	\item 在midstrobe点执行误差检测，相比最佳采样滞后
	\begin{figure}[H]
		\centering
		\includegraphics[scale=0.7]{midstrobe later.png}
		\caption{情况三：在midstrobe点误差检测，相比最佳采样滞后}
	\end{figure}
	如上图所示，易知此时的误差$\mu(k)<0$，\textbf{也就表明此时相比完美采样提前了，Gardner环位同步会根据这个误差通过负反馈环路增大采样间隔修正下次采样的时间，迫使$y(k)$回到符号峰值对应的strobe点去。}
	
	\item 在midstrobe点执行误差检测，相比最佳采样提前
	\begin{figure}[H]
		\centering
		\includegraphics[scale=0.7]{midstrobe before.png}
		\caption{情况四：在midstrobe点误差检测，相比最佳采样提前}
	\end{figure}
	如上图所示，易知此时的误差$\mu(k)>0$，\textbf{也就表明此时相比完美采样滞后了，Gardner环位同步会根据这个误差通过负反馈环路减小采样间隔修正下次采样的时间，迫使$y(k)$回到符号峰值对应的strobe点去。}
\end{enumerate}
\textbf{从上面四种情况可以看出，无论是在midstrobe还是在strobe进行误差检测，Gardner环最终都会回归到在strobe点进行误差检测。}

因此，加入了Gardner算法的Farrow插值器的结构如下图所示，Gardner算法通过NCO控制器不断地修正插值时间$m_k,\mu_k$来期望达到最佳采样点。因此接下来将主要介绍NCO控制器的原理以及实现。
\begin{figure}[H]
	\centering
	\includegraphics[scale=0.7]{farrow_structure.png}
	\caption{Gardner环结构}
\end{figure}
数控振荡器(NCO)在这里起着一个相位递减器的作用，它的差分方程可以表示为，相当于对$[\eta(m)-\omega(m)]$取绝对值
\begin{equation}
	\eta(m+1)=[\eta(m)-\omega(m)]mod1
\end{equation}
式中，$\eta(m)$为第$m$个工作时钟的NCO寄存器的内容，$\omega(m)$为第$m$个工作时钟的控制字，两者都是正小数。

需要注意的是，NCO的工作周期为$T_s$，而内插器的工作周期为$T_i$，$\omega(m)$由环路滤波器进行调节，以使NCO能够在最佳采样点刚好溢出重置。当环路滤波器稳定下来之后，$\omega(m) \approx C_1$，控制字约等于一个常数，所以NCO的值每$T_s$会减小$\omega(m)$，NCO将会每$\frac{T_s}{\omega(m)}$溢出一次。所以$\omega(m)$，内插频率$\frac{1}{T_i}$，采样频率$\frac{1}{T_s}$之间的关系为
\begin{equation}
	\omega(m)  \approx \frac{T_s}{T_i}
\end{equation}
考虑到进入Gardner环的信号采样速率为16MHz，对应每个符号4个采样点，而Gardner环的输出采样速率为8MHz，对应每个符号2个采样点，所以稳定下来后$\frac{T_s}{T_i}=0.5$。这样子就控制到了基本上每2个$T_s$就会让NCO溢出内插一次。

NCO寄存器的内容随时间变化如下图所示
\begin{figure}
	\centering
	\includegraphics[scale=0.7]{Gardner nco.png}
	\caption{NCO寄存器内容变化图}
\end{figure}
图中$m_kT_s$为第$m_k$个时钟脉冲，在$(m_k+1)T_s$个时钟脉冲的时候，NCO的值为$\eta(m)-\omega(m)$，发现这个时候NCO的值小于0，也就是发生了一次溢出，因此和1进行$mod$运算的时候会变成当前NCO的绝对值。NCO发生了一次溢出意味需要在这个$m_kT_s$和$m_kT_s$之间内插一个点进去，根据相似三角形可以计算出
\begin{equation}
	\frac{\mu_k}{\eta(m_k)} = \frac{(1-\mu_k)T_s}{1-\eta(m_k+1)}
\end{equation}
解出分数间隔$\mu_k$为
\begin{equation}
	\mu_k = \frac{\eta(m_k)}{\omega(m_k)}
\end{equation}

所以整个Gardner环位同步的MATLAB迭代流程为：
\begin{enumerate}
	\item 获取最新的采样时刻$(m+1)T_s$的数据，计算当前NCO寄存器的值$\eta(m+1)=\eta(m)-\omega(m)$，观察是否发生了溢出，如果发生了溢出需要取NCO寄存器的绝对值并计算分数间隔$\mu_k$，并且标注此处发生了溢出；
	\item 根据分数间隔$\mu_k$对已有的数据进行Farrow插值
	\item 根据Gardner误差检测算法，需要两个strobe点才能计算误差，也就是每插值两次计算一次定时误差
	\item 根据定时误差经过环路滤波器修正NCO的控制字$\omega$，从而间接控制分数间隔$\mu_k$
\end{enumerate}

完成Gardner环位同步之后，由于每个符号有两个采样点(strobe和midstrobe)，而我们只需要strobe，因此我们只需要按照每两个Gardner环位同步输出取一个就可以取得最后的符号采样值进入下一个环节。




\subsection{均衡}\label{subsec:software equalization}
考虑到多径信道的影响，信号经过多径信道传输之后在时域和频域上都会有显著的变化。如果不经过均衡直接对信号进行解调等操作，其结果如下图所示，可以看到此时的星座图是十分杂散的，毫无规律的，这个时候解调出来的数据显然不是发送时的数据，也不是我们期望的数据。因此我们需要对接收到的信号进行均衡，也就是将这个杂散的星座图还原成QPSK标准的星座图。

{\color{red} 缺图}


因此，为了消除非理想信道传输带来的影响，需要对信号进行均衡。\textbf{首先需要明确的是均衡处于QPSK收发机的什么位置，无论是基于导频训练的均衡还是盲均衡都是作用在符号上的，因此均衡需要输入QPSK接收机接收到的符号，考虑到一个完整的QPSK收发机的组成包括信源-采样-成型-调制-信道-解调-下采样-Gardner最佳采样-PLL等，因此均衡应该放在Gardner环最佳采样之后。}

均衡部分结构图如下所示，可以看到均衡滤波器的输入是经过信道传输之后的\textbf{符号}，而期望的则是发送的\textbf{符号}，一般的均衡滤波器使用的都是梯度下降算法，令滤波器的系数$\textbf{W}$不断迭代使滤波器的输出\textbf{符号}和发送的\textbf{符号}之间的误差尽可能地小。
\begin{figure}[H]
	\centering
	\includegraphics[scale=0.4]{equalizer_structure.png}
	\caption{均衡模块系统框图}
\end{figure}

针对是否需要导频作为参考可以分为基于导频训练的CLMS均衡算法以及不需要导频训练的CMA盲均衡，接下来将分开介绍两种算法。
\subsubsection{基于导频训练的CLMS均衡}
对于实数均衡LMS(Lease-Mean-Squre)算法，$j$时刻的输入矢量$X_j$以及均衡滤波器系数$W_j$可以表示为：
\begin{equation}
	\begin{split}
		X_j &= \begin{bmatrix}
			x_1,x_2,x_3,\cdots,x_j
		\end{bmatrix}^T\\
		W_j &= \begin{bmatrix}
			w_1,w_2,w_3,\cdots,w_j
		\end{bmatrix}^T
	\end{split}
\end{equation}
$j$时刻的输出可以表示为
\begin{equation}
	y_j = W_j^TW_j = W_j^TX_j
\end{equation}
$j$时刻的误差则被定义为期望输出$d_j$以及滤波器输出$y_j$之间的差值
\begin{equation}\label{eq:error}
	\epsilon_j = d_j-y_j=d_j-W_j^TX_j
\end{equation}

而对于复数信号处理而言，由于输入信号$x_j$为复数，LMS算法需要适应复数域的计算，因此在LMS的基础上修改为CLMS(Complex LMS)算法。

定义输入信号矢量$X_j$以及滤波器系数矢量$W_j$为
\begin{equation}
	\begin{split}
		X_j&=
		\begin{Bmatrix}
			x_{1Rj}\\
			x_{2Rj}\\
			\vdots\\
			x_{nRj}	
		\end{Bmatrix}
		+i\begin{Bmatrix}
			x_{1Ij}\\
			x_{2Ij}\\
			\vdots\\
			x_{nIj}	
		\end{Bmatrix}
		=X_{Rj}+iX_{Ij}\\
		W_j&=
		\begin{Bmatrix}
			w_{1Rj}\\
			w_{2Rj}\\
			\vdots\\
			w_{nRj}	
		\end{Bmatrix}
		+i\begin{Bmatrix}
			w_{1Ij}\\
			w_{2Ij}\\
			\vdots\\
			w_{nIj}	
		\end{Bmatrix}
		=W_{Rj}+jW_{Ij}\\
	\end{split}
\end{equation}
其中$R$表示实数部分，$I$表示虚数部分。复数信号误差$\epsilon_j$以及期望信号$d_j$都得以复数的形式表示：
\begin{equation}
	\begin{split}
		\epsilon_j&=\epsilon_{Rj}+i\epsilon_{Ij}\\
		d_j&=d_{Rj}+id_{Ij}
	\end{split}
\end{equation}
均衡滤波器的输出$y_j$同样为一复信号，表示为：
\begin{equation}
	y_j=y_{Rj}+iy_{Ij}
\end{equation}

CLMS算法必须能够同时修改均衡滤波器系数$W_j$的实部和虚部部分，在某种程度上使误差的实部和虚部同时最小。需要注意的是，考虑到实部和虚部为正交关系，将他们独立分别最小化是不可行的。一个合理的目标是使平均总误差功率最小，如下式所示：
\begin{equation}
	E[\epsilon_j\bar{\epsilon_j}]=E[\epsilon_{Rj}^2+\epsilon_{Ij}^2]=E[\epsilon_{Rj}^2]+E[\epsilon_Ij^2]
\end{equation}
其中$E$表示数学期望，$\bar{\epsilon}$表示误差矢量$\epsilon$的共轭。

CLMS算法的推导过程和LMS算法的推导过程类似，只是需要考虑复数域的运算。式\ref{eq:error}的共轭表示为：
\begin{equation}
	\bar{\epsilon_j}=\bar{d_j}-\bar{W_j}^T\bar{X_j}=\bar{d_j}-\bar{X_j}^T\bar{W_j}
\end{equation}
平均总误差功率关于滤波器系数实部的梯度可以表示为：
\begin{equation}
	\nabla R(\epsilon_j\bar{\epsilon_j})=
	\begin{Bmatrix}
		\frac{\partial(\epsilon_j\bar{\epsilon_j})}{\partial W_{1R}}\\
		\vdots\\
		\frac{\partial(\epsilon_j\bar{\epsilon_j})}{\partial W_{nR}}
	\end{Bmatrix}
	=
	\epsilon_j\nabla R(\bar\epsilon_j)+\bar{\epsilon_j}\nabla R(\epsilon_j)
	=\epsilon_j(-\bar{X_j})+\bar{\epsilon_j}(-X_j)
\end{equation}
同理，平均总误差功率关于滤波器系数虚部的梯度可以表示为：
\begin{equation}\label{eq:ewr}
	\nabla R(\epsilon_j\bar{\epsilon_j})=
	\epsilon_j\nabla I(\bar\epsilon_j)+\bar{\epsilon_j}\nabla I(\epsilon_j)
	=\epsilon_j(-\bar{iX_j})+\bar{\epsilon_j}(-iX_j)
\end{equation}
更新滤波器系数的实部和虚部，则有
\begin{equation}\label{eq:ewi}
	\begin{split}
		W_{Rj+1}=W_{Rj}+\mu\nabla R(\epsilon_j\bar{\epsilon_j})\\
		W_{Ij+1}=W_{Ij}+\mu\nabla I(\epsilon_j\bar{\epsilon_j})
	\end{split}
\end{equation}
考虑到滤波器系数$W_j=W_{Rj}+iW_{Ij}$，则复数形式的迭代公式可以写成：
\begin{equation}
	W_{j+1}=W_j-\mu[\nabla R(\epsilon_j\bar{\epsilon_j})+i\nabla I(\epsilon_j\bar{\epsilon_j})]
\end{equation}
联立式\ref{eq:ewr}和式\ref{eq:ewi}有：
\begin{equation}
	W_{j+1} = W_j+2\mu\epsilon_j\bar{X_j}
\end{equation}

因此MATLAB的迭代过程可以如下表示:
\begin{enumerate}
	\item 读取$j$时刻的接收机接收信号，求共轭，存入输入信号矢量$\bar{X_j}$中
	\item 对输入信号矢量滤波求得均衡滤波器输出$y_j=W_j^TX_j$
	\item 求期望信号$d_j$和滤波器输出$y_j$之间的误差$\epsilon_j=d_j-y_j$
	\item 迭代更新滤波器的权重$W_{j+1} = W_j + 2\mu\epsilon_j\bar{X_j}$
\end{enumerate}

需要注意的是：为了实现上述迭代，在第3步中需要期望信号$d_j$来计算误差$\epsilon_j$，这个期望信号的集合被称为训练序列，也就是说发射机需要发送训练序列来供接收机更新均衡滤波器的系数。

训练序列不含任何信息，占用了大量的带宽，因此降低了通信系统的效率。除此以外，由于之前环节的滤波器、信道延时等带来的符号延时，如何确定当前接收到的符号对应训练序列的哪一位需要使用帧同步环节。


\subsubsection{CMA盲均衡}
为了减少通信带宽的消耗，一种不依赖训练序列的盲均衡方法被提出，这种仅用接收信号本身的先验信息来均衡信道特性，使其输出序列尽量逼近发送序列的自适应均衡技术，能够有效地补偿信道的非理想特性，克服码间串扰，减小误码率，提高通信质量。

\textbf{对于频率调制和相位调制信号，如QPSK信号，其包络为常数}，利用该性质，可以使均衡器的输出的模值与规定的某个常数值之间的距离最小化，以达到均衡的目的。\textbf{其核心思想是利用信号的模值统计信息，代替发送信号的训练序列来构造代价函数，通过最小化代价函数推导出均衡器抽头系数的更新公式，不断迭代，使均衡器达到最优状态。}这里假设接收机接收到的符号为$x$，均衡滤波器的系数为$h$，均衡滤波器的输出为$y$，则可以定义代价函数$J$为：
\begin{equation}
	\begin{split}
		J &= E[e_{Godard}^p(k)]\\
		&= E[(|y(k)|^q-r_q)^p]\\
		&= E[(|w^H(k)x(k)|^q-r_q)^p]
	\end{split}
\end{equation}
其中$r_q=\frac{E[|x(k)|^{2q}}{E[|x(k)|^q]}$，$p$和$q$都是正整数。$r_q$的值定义了$|y(k)|^q$应该达到的大小，同时对它施加一个幂为$p$的惩罚。与上面用到的CLMS算法相似，这里也用到了瞬时梯度作为估计值，也就是对上面式子中的$w^*(k)$求微分如下：
\begin{equation}
	\begin{split}
		\nabla_{w^*}J &= pe_{Godard}^{p-1}(k)\frac{\partial (|y(k)|^q)}{\partial w^*(k)}\\
		&= \frac{1}{2}pe_{Godard}^{p-1}(k)(|y(k)|^{q-2})\frac{\partial (|y(k)|^2)}{\partial w^*(k)}\\
		&= \frac{1}{2}pe_{Godard}^{p-1}(k)(|y(k)|^{q-2})\frac{\partial (w^*(k)x(k)x^H(k)w(k))}{\partial w^*(k)}\\
		&= \frac{1}{2}pe_{Godard}^{p-1}(k)(|y(k)|^{q-2})x(k)X^H(k)w(k)\\
		&= \frac{1}{2}pe_{Godard}^{p-1}(k)(|y(k)|^{q-2})y^*(k)x(k)\\
	\end{split}
\end{equation}
这里涉及到的相关的矩阵运算公式有：
\begin{equation*}
	\begin{split}
		y(k)&=w^H(k)x(k)\\
		y^*(k)&=w^T(k)x^*(k)=x^H(k)w(k)\\
		\nabla_{z^*}(Z^HRZ)&=RZ\\
	\end{split}
\end{equation*}
可以得到随机梯度形式的Godard算法，其迭代方程为：
\begin{equation}\label{eq:Gadard iter}
	w(k+1) = w(k)-\frac{1}{2}\mu pq(|y(k)|^q-r_q)^{p-1}|y(k)|^{q-2}y^*(k)x(k)
\end{equation}
当$p=q=2$时，就是常模CMA算法，MATLAB的迭代过程为：
\begin{enumerate}
	\item 初始化$x(0)$和$w(0)$为随机向量
	\item $r_2=\frac{E[|x(k)|^{4}]}{E[|x(k)|^2]}$，通常对于标准的QPSK而言，$r_2=1$
	\item 滤波处理$y(k)=w^H(k)x(k)$
	\item 计算误差$err_{CMA}(k)=|y(k)|^2-r_2$
	\item 更新滤波器系数$w(k+1)=w(k)+2\mu err_{CMA}(k)y^*(k)x(k)$，\textbf{需要注意的是，和}\ref{eq:Gadard iter}\textbf{相比，这里的$err_{CMA}(k)=y(k)(|y(k)|^2-r_2)$}
\end{enumerate}

CMA的算法根本目的是使信号的模尽可能地等于1，也就是使信号尽可能地收敛到单位圆上，但是CMA只能做到这种程度了，收敛到单位圆上并不能确保相位一定满足QPSK调制时的相位，这个时候存在着一个相偏的问题，需要加入帧同步以及PLL锁相环环节来解决相偏的问题。
\subsubsection{均衡结果}
完成均衡后的星座图如下图所示，可以看到相比于均衡前的结果，均衡后的结果呈现出来典型的QPSK星座图，虽然相位没有对齐，这也是非理想信道传输对信号带来的相偏频偏的影响。

{\color{red} 缺图}




\section{硬件设计}
硬件设计总体上分为发射机设计和接收机设计，考虑到不同模块的时钟不同，发射机的硬件设计如下所示：
\begin{figure}[H]
	\centering
	\includegraphics[width=1\textwidth]{fpga_sender_structure.png}
	\caption{发射机硬件设计结构图}
\end{figure}

信源是向发射机提供需要发送的数据的设备，考虑到信源发送数据存在时域上的不确定性，每次发送信号的长度、间隔时长都不确定，并且发送给发射机的速率通常大于发射机设计的符号速率，因此需要用一个缓冲区Buffer(FIFO)完成数据的缓存，发射机从缓冲区中不断的读取数据完成后续的处理。

数据一旦从缓冲区中取出将会\textbf{连续不断地}经过若干个模块的处理，和软件设计发射机部分相同，发射机首先以5MHz的时钟从缓冲区中不断地取出数据形成比特流，QPSK符号映射需要两个比特映射为一个符号，因此符号映射后形成的符号速率(码率)为2.5 symbols/s，IQ通道每个符号会被4倍过采样，因此上采样部分需要10MHz的时钟驱动。经过根升余弦滤波器成型滤波之后的IQ通道数据会被再次6倍过采样形成发送给DAC的数据，因此过采样部分需要60MHz的时钟驱动。

\textbf{考虑到FPGA最后向AD9361传输的是12bits的数据，因此在这里设计硬件实现中的数据以12bits为标准，后续需要运算的时候对12bits进行对应的符号位扩展。}

\subsection{QPSK映射}\label{subsec:hardware qpsk mapping}
根据QPSK符号映射软件部分\ref{subsec:software qpsk mapping}的工作可知，硬件实现QPSK映射的关键在于需要读取两个bits对他们进行按照规律的映射。可以采用计数器实现，当计数器从0开始计数重新计数到0时就可以进行一次映射。

\subsection{升余弦成型滤波器}\label{subsec: hardware rcos shapping}
根据软件部分\ref{subsec:software rcos shapping}的工作，对于升余弦滤波器的硬件实现，分为两个部分，一个是对信号进行4倍过采样，另一个部分是对信号进行升余弦滤波器滤波处理。

其中4倍过采样参照MATLAB中的采样方式，在每两个数据点中插入3个0，因此可以使用计数器实现。在4倍数据时钟驱动下，令计数器从1开始计数，只有当计数器计数到4时才用数据更新输出的数值，其他时刻都输出0。

而升余弦滤波器可以通过MATLAB的rcosdesign函数预先设计出来，然后对设计出来的升余弦滤波器的系数进行定点化处理，可以采用MATLAB的filterDesigner提供的定点化功能实现，也可以手动计算实现。将定点化之后的滤波器系数导出为.coe文件，在vivado中新建一个blockdesign，导入vivado提供的fir的ip核，设置ip核的Filter Coefficients中的Coefficients File为MATLAB生成的.coe文件，可以在左侧的Freq Response窗口看到这个滤波器的频域响应。

\textbf{需要注意的是，由于MATLAB生成的滤波器系数都是小数，而FPGA需要的滤波器系数都是整数，因此定点化之后的滤波器相比MATLAB原始的滤波器有了定点化带来的增益。以16bits的有符号数定点化为例，定点化之后相比定点化之前多了$20log_{10}{2^{15}}=90.3096dB$的增益。}

设置FIR滤波器的输入为12bits，而FIR中的系数都是16bits的定点化小数。同时在Channel Specification中设置对应的时钟Clock Frequency以及Input Sampling Frequency为驱动FIR运转的时钟，此处为10MHz。最终vivado中的FIR滤波器的IP核设置如下图所示：
\begin{figure}[H]
	\centering
	\includegraphics[width=1.0\textwidth]{fpga_rcos_setting.png}
	\caption{vivado中的FIR滤波器的IP核设置}
\end{figure}

4倍采样以及升余弦滤波器的modelsim仿真如下图所示，以i通道的数据为例，由上到下依次为输入信号data\_i，4倍过采样之后的信号data\_4x\_i，升余弦成型的结果temp\_shaping\_i，以及最后截位后输出的结果data\_rcos\_i。可以观察到4倍过采样在两个数据之间插入0使得采样后的波形相比采样前的波形更细，成型之后的波形符合MATLAB中仿真的结果。
\begin{figure}[H]
	\centering
	\includegraphics[width=1\textwidth]{fpga_sample_4x_and_rcos_shaping.png}
	\caption{升余弦modelsim仿真结果}
\end{figure}

\textbf{这里需要注意的是，成型之后的结果temp\_shaping\_i是FIR的输出，是31bits的数据，而我们期望的最终输出是12bits的，因此这里需要进行截断。考虑到输入的12bits中都是整数位，FIR的参数中15bits为小数位，因此temp\_shaping\_i中有15位小数位，为了使符号位长度一致，这里截断最后4位小数，输出11个小数位。}


\subsection{上采样}\label{subsec:hardware upsample}
参照软件设计中\ref{subsec:software upsample}的工作可以知道，要实现Farrow内插器，本质上就是实现式\ref{eq:farrow}的公式，如下所示：
\begin{equation}\label{eq:farrow}
	\begin{split}
		f_1 &= 0.5x(m)-0.5x(m-1)-0.5x(m-2)+0.5x(m-3)\\
		f_2 &= 1.5x(m-1)-0.5x(m)-0.5x(m-2)-0.5x(m-3)\\
		f_3 &= x(m-2)\\
		y_1(k) &= f_1 u(k) u(k) + f_2 u(k) + f_3\\
	\end{split}
\end{equation}

从公式中可以看出，要实现一个Farrow内插器，需要以下元素，首先是4个输入信号$x(m), x(m-1), x(m-2), x(m-3)$，不断变化的内插的分数间隔$\mu(k)$。\textbf{需要注意的是，在MATLAB中不需要考虑时序问题，而在FPGA实现的时候，加法以及乘法运算都是需要时间的，因此加法以及乘法完成之后的当前的分数间隔$\mu(k)$已经不是当时运算之前的分数间隔了，必须要对相应的变量进行延时寄存。这样在进行下一步运算的时候才可以用正确的时序对应的变量进行计算。}

这里设定加法器和乘法器都为pipeline stage为1的结构，也就是普通的时序逻辑，在下一个上升沿的时候会输出运算的结果，因此这里需要注意\textbf{产生了一个延时。}整个时序可以如下分解：
\begin{enumerate}
	\item 在第一个上升沿，同时计算$0.5x(m)-0.5x(m-1)$与$0.5x(m-2)-0.5x(m-3)$，$1.5x(m-1)$与$0.5x(m)+0.5x(m-2)$。对于$1.5x(m-1)$可以使用移位运算实现；
	\item 在第二个上升沿，加法器已经完成了上述的所有运算，因此可以进行$f_1$和$f_2$的运算；
	\item 在第三个上升沿，$f_1$和$f_2$都已经完成计算，可以进行乘法运算，这里使用两个乘法器同时计算$f_1*u(k)$和$f_2*u(k)$；
	\item 在第四个上升沿乘法器完成了上述的计算，这里使用一个乘法器计算$f_1*u(k)*u(k)$，由于需要在进行一次乘法运算，因此需要对$f_2*u(k)$进行一次延时，以确保在下一个上升沿的时候需要计算的元素是对齐的。
	\item 在第五个上升沿的时候，按道理需要完成最后的运算$y_1(k) = f_1 u(k) u(k) + f_2 u(k) + x(m-2)$。但是这里需要考虑到\textbf{最后一步计算中有$x(m-2)$这个变量}，当前上升沿的$x(m-2)$和五个上升沿前的$x(m-2)$有什么区别是值得分析的问题。这里有两种情况，第一种是5个上升沿前的信号$x(m-2)$和当前的信号相等，比如一个信号的第1个内插点和第6个内插点；另一种是不相等的情况，比如一个信号的第4个内插点和下一个信号的第3个采样点。如果简单的认为只有一种情况势必会导致另一种情况的确实，从而导致信号的不连续，因此这里需要对信号进行统一的延时确保只有第二种情况会发生。由于此处的Farrow内插倍数为6，因此这里对$f_1 u(k) u(k)$和$f_2 u(k)$都进行2个clk的延时，这样整体的时延就变成了7个clk，确保当前时刻应该用$x(m-3)$来参与计算。
	\item 在第七个上升沿的时候，输出整体计算结果。
\end{enumerate}

明确了思路了之后，在编写verilog代码的时候需要注意到定点数运算的一些细节问题。因为farrow内插牵扯到加法和乘法运算，因此需要注意有效数据位的问题。在\ref{subsec: hardware rcos shapping}中，我们规定了最后的输出为12bits，由包含1个表示符号的整数位以及11个表示小数的bit，为了防止后续运算可能导致的溢出问题，这里对数据进行符号位扩充到16bits，也就是5个整数位，11个小数位。farrow内插时的乘法运算的小数uk设置为16bits，1个符号位+15个小数位。因此进行乘法运算的结果为32bits，包括6个整数位以及26个小数位。但是为了进行接下来的乘法运算，需要依旧保留11个小数，因此需要对乘法运算的结果截断为11个小数，所以舍弃乘法结果的低15的bits。

farrow内插器整体的增益为1，设计最后输出到DAC为12bits，因此最后需要对扩充到16bits的数据进行最终的截断，截断最高4个bits作为输出给DAC的数据。

farrow内插器的modelsim仿真结果如下图所示，可以看到相比于内插前的波形，farrow内插之后的波形更加的顺滑，并且没有突变等情况。
\begin{figure}[H]
	\centering
	\includegraphics[width=1.0\textwidth]{fpga_farrow_comparision.png}
	\caption{farrow内插器modelsim仿真结果}
\end{figure}

\section{仿真及实验结果}

\section{总结}

\section{问题汇总}
本部分将记录在进行程序编写仿真时遇到的问题以及搜索的解决办法以及相关资料。
\begin{enumerate}
    \item \textbf{Vivado Error: (vlog-7) Failed to open design unit file.}
        
    一般出现此类问题是因为Tcl console中提到的对应的某个block的名字过长导致无法读取该文件，尝试重新创建一个block取一个短一点的名字。
    
    \href{https://programmerah.com/vivado-error-vlog-7-failed-to-open-design-unit-file-43918/}{Vivado Error: (vlog-7) Failed to open design unit file}
    
\end{enumerate}



\bibliographystyle{IEEEtran}
\bibliography{ref} 
\end{document}